Entanglement entropy and entanglement spectrum of the Kitaev model 
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In this paper, we obtain an exact formula for the entanglement entropy of the ground state and 
all excited states of the Kitaev model. Remarkably, the entanglement entropy can be expressed 
in a simple separable form S = So + Sf, with Sf the entanglement entropy of a free Majorana 
fermion system and So that of a Z2 gauge field. The Z2 gauge field part contributes to the universal 
"topological entanglement entropy" of the ground state while the fermion part is responsible for the 
non-local entanglement carried by the Z2 vortices (visons) in the non-Abelian phase. Our result also 
enables the calculation of the entire entanglement spectrum and the more general Renyi entropy 
of the Kitaev model. Based on our results we propose a new quantity to characterize topologically 
ordered states — the capacity of entanglement, which can distinguish the states with and without 
topologically protected gapless entanglement spectrum. 



Exotic phases such as fractional quantum Hall (FQH) 
states, which are not in the paradigm of conventional 
symmetry breaking, were termed as topologically ordered 
[l| since they have robust ground state degeneracy which 
is protected against all local perturbations, but sensitive 
to the topology of the system A topologically or- 
dered state has non-local pattern of quantum entangle- 
ment, which is essential for the proposal of topological 
quantum computation 0-0]. 

By bipartitioning a system spatially, the entanglement 
entropy (EE) measures how closely entangled the two 
subsystems are. For a gapped system, EE is usually pro- 
portional to the area of the interface between the two 
subsystems in the thermodynamic limit. However, as 
discovered ^, ^ by Levin and Wen as well as Kitaev and 
Preskill, the entanglement entropy of a topologically or- 
dered state contains a universal constant term, which is 
uniquely determined by the topological order of the state, 
named as topological entanglement entropy (TEE) . TEE 
enables a direct characterization of topological ordered 
states without referring to the Hamiltonian. EE and TEE 
are properties of a many-body state and are usually hard 
to compute. EE and/or TEE have been computed ex- 
actly or numerically for several models such as toric code 
model^ayShja], FQH states[lM2| and quantum dimer 
models [13Hl5| . So far there has been no exact result for 
the EE of topologically ordered states whose quasiparti- 
cles obey non-Abelian statistics. 

This paper serves to fill in that gap by providing a 
simple but exact method to compute the EE for any 
eigenstate (either ground or excited states) of the Ki- 
taev model [l6| . which is one of the most important 
exact solvable models with non-Abelian anyons. The 
essence of our method is a rigorous proof that the EE 
of the Kitaev model is equal to that of two decoupled 
systems: a sourceless Z2 gauge field and a free Majorana 
fermion system. Although the TEE of the ground state 



comes only from the Z2 gauge field, the fermionic part 
is responsible for all nontrivial entanglement properties 
of the non-Abelian phase. Besides the EE, our method 
also enables the computation of the whole entanglement 
spectrum (ES), i.e., the eigenvalue spectrum of the re- 
duced density matrix 17[. We show that the entangle- 



ment spectrum is gapless or gapped in the non-Abelian 
and Abelian phase of the Kitaev model, respectively. We 
propose a new quantity, the capacity of entanglement, 
which can be used to distinguish different topological 
states with gapped and gapless entanglement spectrum. 

Kitaev model is a spin- 1/2 model originally proposed 
on the honeycomb lattice [16| with the Hamiltonian 
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where x, y, z-link stand for the three types of links. It has 
a non-Abelian phase when the time-reversal symmetry is 
broken either explicitly by magnetic field jl6l | or three- 
spin couplings [18], or spontaneously by decorating the 
honeycomb lattice [l^ . For simplicity, hereafter we will 
present our exact results of EE and entanglement spec- 
trum for the Kitaev model on honeycomb lattice, but our 
approach can be generalized straightforwardly to a broad 
class of Hamiltonians, including the Kitaev model on any 
trivalent lattice and Gamma matrix models (iol - fiij . 

The Kitaev model can be solved by introducing the 
Majorana representation of the Pauli matrices [IB]: erf = 
*7i"^i = z), where 7" and 77^ are Majorana 
fermion operators. 7^ and rji on each lattice site de- 
fine a 4-dimensional Hilbert space, so that the Majo- 
rana representation of a spin-1/2 is redundant. The 
physical Hilbert space is defined by a constraint Di = 
—iafa^af = "yfl^liTji — 1. In other words, a state \'^) is 
physical only if Di\^) = \^) for every i. In the Majorana 
representation we have erf cr" — jfjfrjirjj — —iuijrjirjj, 
in which the link operators uij — i^f^f mutually com- 
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mute and also commute with the Hamiltonian. Since 
ufj = 1, Uij can be considered as c- numbers with values 
Uij — ±1, so that the Kitaev model is equivalent to a 
free model o f 77 Ma.jorana fermions coupled to static Z2 
gauge fields [l6l.[llfl9l.[23-[27t. The ground state of such 
a model is given by the direct product of a Z2 gauge con- 
figuration \u) and the corresponding Majorana fermion 
ground state \4>{u)). Here the configuration u is deter- 
mined by minimizing the fermion ground state energy. 
There is a macroscopic ground state degeneracy in the 
enlarged Hilbert space, because each state \u) \4>{u)) is 
degenerate with all the states \u') (g) \(j>{u')) with u' gauge 
equivalent to u. However, such a degeneracy is removed 
when the constraint = 1 is applied. The physical 
ground state is the "gauge" average of the degenerate 
states, implemented by the projection 16|: 
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where N is the total number of lattice sites, g denotes a 
set of lattice sites, and Dg = Yl^^gDi- We define D = 
riiGC ■^ith £ the set of all lattice sites. The sum 
is taken over all possible subsets g oi £. Note that, in Eq. 

we implicitly assumed that D \u) (g) = \u) <8) 

\ip{u)) because states with D = —1 will be annihilated 
by the projection. Consequently, we have Dg = DDg for 
the complement g — £ — g, so that Dg \u) (g) \(l>{u)) = 
Dg=c-g \u) ^ In other words there are only 2-^^^ 

inequivalent gauge transformations, as expected. 

We define the Kitaev model on a torus and bipar- 
tite the lattice into subsystems A and B, as shown 
in Fig. [TJ The EE between A and B is defined as 
S = -TrA [pAlogpA], where pA = Tr^p = Tr^ 1^-) {^\ is 
the reduced density matrix of A. To calculate the EE, 
we will follow the replica tridi^ introduced in Ref. [2^ 



S = -Tta [pa log pa] = "I^TrA [pV\ 



(3) 



The entanglement entropy can be obtained if we can com- 
pute Tr^ [pa] arbitrary positive integer n and then 
extrapolate the result to n g R. 




FIG. 1: The schematic honeycomb lattice is bipartitioned into 
two parts A and B. The partition boundary (dashed line) cuts 
the links Unbn, n — 1, • ■ • , 2L. New Z2 gauge variables (see 
text) WA.n and WB,n are introduced on the new (dotted) links 
a2n-ia2n and 6271-162™, n = 1, ■ ■ ■ , L, respectively. 



To obtain pA, we trace out the spin degree of free- 
dom in B, which normally can be carried out in terms 
of fermions and gauge fields. However, the gauge fields 
on the partition boundary are shared by A and B; so 
we regroup those gauge fields on the boundary links to 
introduce new Z2 gauge variables which lives in A and B 
exclusively, as shown in Fig. [TJ (see supplementary ma- 
terial [i^l for details.) The calculation of Tr[pA] requires 
some careful treatment of the gauge transformation but 
is a well-defined mathematical procedure. Thus, we will 
leave the details involved in obtaining pA and Trp^ to 
the supplementary material [29] and present only the fi- 
nal result here: 



TrA[pl]=TrA,G[pl,G]-TrA,F[pV], 



(4) 



for any positive integern. Here pa,f = TrB[|(?!>(u)) 
and pa.g — T^^b[\G{u)) {G{u)\] are the reduced density 
matrices for the free Majorana fermion state \(j){u)), and 
a pure Z2 gauge field [9|, respectively, and the ground 
state of the Z2 gauge field \G{u)) is given by a equal 
weight superposition of all the 2-^^"^ gauge field con- 
figurations \u) that are gauge equivalent to \u), i.e., 
\G{u)) = 2-(^-i)/2X;i^„|u). Physically, such a sim- 
plification occurs because the effect of the gauge trans- 
formation Dg on the fermion state in region B is canceled 
out once the trace over gauge field configurations is taken. 

Combining Eq. ([4]) and Eq. ([3]) , it is now obvious that 
the EE S can be separated into gauge field part Sg and 
fermion part Sp as follows: 



S 



(5) 



Eq. ([4]) and jS]) are among the central results of this work. 
By explicit calculation j29| one can obtain Tr^, g[Pa g\ ~ 
2-(^-i)("-i),sothatS'G = (L-l)log2. As wiU be shown 
below, the fermion part has the form Sp = otL + o(l), 
where a is a positive constant and oil) represents terms 
which vanish as L — >■ 00. In the thermodynamic limit, 
the total entanglement entropy is given by 



5*= (a -flog 2)L- log 2, 



(6) 



from which we conclude that the TEE is S'topo = — log 2. 
Our derivation is valid for all phases of Kitaev model, 
including the Abelian {Z2 gauge theory), non-Abelian 
(Ising anyon) phases, and also gapless phases. Thus our 
result directly proves that the TEE for the Abelian and 
non-Abelian phases are identical, as expected from the 
total quantum dimensions of their quasiparticles (30l | . 

Despite its trivial contribution to TEE, the fermion 
sector Sp is responsible for all the essential differences 
between the Abelian and non-Abelian phase of Kitaev 
model in their quantum entanglement properties. The 
EE of a free fermion system can be computed by the 
method introduced in Ref. 3l|. To obtain an explicit 



understanding to the fermion EE we consider a torus 
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divided by two parallel circles into A and B regions, as 
shown in Fig. [2] (a). The boundary circle is along the 
y direction. On the torus, the free Majorana fermion 
Hamiltonian can be block-diagonalized in the basis of 
ky-. H = Y.r.jmV3hi3 = J2x,x'.k^vUky)'T-^^'i>^y)Vx'{ky). 
Thus the system consists of decoupled one-dimensional 
subsystems of each ky . The EE is given by [sij 

sf^-Ij2 log + (1 - iog(i - ^»)] (^y)' 

where \n{ky) are the eigenvalues of the single-particle 
correlation function Cxx'{ky) = {'ill{ky)r]x'{ky)) for 
each ky. A„ plays the role of Fermi-Dirac distribution 
l/(e^'^" _j_ 2) in thermal entropy, so that A,i = (1) cor- 
responds to fully occupied (unoccupied) states, respec- 
tively. The "entanglement spectrum" (ES) Xn{ky) has 
been computed numerically for both non-Abelian and 
Abelian phases, as shown in Fi g, p (b), for the Kitaev 
model with three-spin terms J' [ISj. The ES is gapped 
for the Abelian phase, and gapless for the non-Abelian 
phase, similar to the edge states in the energy spectrum. 
Similar observation has been made in topological insu- 
lators and superconductors 32-3^ and in FQH systems 
17| . The two gapless branches in the ES come from 



the two boundaries between A and B. Since A„(A;y)'s 




FIG. 2: (Color online) (a) Schematic picture of a torus and 
a cylinder, each split to two regions A and B. The cylinder 
is equivalent to a sphere with two quasi-particles. (b) The 
entanglement spectrum Xn{ky) versus ky for non-Abelian (red 
solid lines) and Abelian (blue dotted lines) state on torus. 
Here and below, we take the parameters Jx ~ Jy ~ Jz ~ 
1 and next-nearest neighbor coupling J' — 0.2 for the non- 
Abelian state, and Jx = Jz = 1, Jy ~ 2.5, J' = 0.2 for the 
Abelian state, (c) The entanglement spectrum for the non- 
Abelian state on cylinder. The blue circle marks an additional 
state with A — 1/2 at ky = 0. (d) The entropy S{ky) versus 
ky for non-Abelian (red solid line) and Abelian (blue dotted 
line) states on the torus and for non-Abelian state on cylinder 
(black dashed line with circles). 



are smooth functions of ky, we see from Eq. ([7]) that 
in the continuum limit Sp — S{ky) ~ L J S{ky)^^ 
satisfies the area law. It is interesting to note that a 
"gap" always exists between the edge states and other 
bulk states with A„(fcj,) close to or 1, which is analo- 
gous to the "entanglement gap" studied in Ref. [s^l for 
FQH system. 

The situation becomes more interesting when we con- 
sider a cylinder with periodic boundary condition (PBC) 
and the partition shown in Fig. [2] (a). As shown in 
Fig. [21(c), in the non-Abelian phase the numerical cal- 
culation gives only one branch of "gapless" states in the 
entanglement spectrum. Physically, this is because the 
coupling through the other boundary between A and B 
is removed by the open boundary condition. However, at 
ky = there is one isolated additional state with A = 1/2, 
as shown by the blue circle in Fig. [2] (c), which is due 
to the non-local entanglement between the two Majorana 
zero modes at the open boundary. Consequently, the en- 
tropy S{ky) is not a smooth function of ky but has an 
additional log \/2 contributed by ky — 0, as shown in 
Fig. [21(d). Compared with the torus case, in the ther- 
modynamic limit we get Sp — aL + log -^72, which shows 
explicitly that in the non-Abelian phase a cylinder with 
PBC is topologically equivalent to a sphere with two non- 
Abelian quasiparticles (usually named as a particles) , as 
illustrated in Fig. [21 (a). Each particle carries a log a/2 
entropy which is solely contributed by the fermion sector. 

Besides the EE, more information is contained in our 
result. The fact that Eq. dH) holds for any posi- 
tive integer n indicates that the many-body entangle- 
ment spectrum — the eigenvalue spectrum of pA is the 
direct product of the ones of pa,g and pa,f- From 
TrA,G [pa^g] = 2-(-^-i)("-i), one'can know 'that pa.g 
has 2^~^ nonzero eigenvalues, all of which are degen- 
erate and have the value of 2"^^"^-'. Consequently all 
non-vanishing eigenvalues of pA are given by those of 
the Majorana fermion reduced density matrix pa,f times 
2~^'^~^\ Thus the low "energy" (i.e., close to the maxi- 
mal eigenvalue of pa) feature in the entanglement spectra 
of Pa can be entirely characterized by its fermionic part, 
which is gapped in the Abelian phase and gapless in the 
non-Abelian phase, as shown in Fig. [21 (b) . 

Such a qualitative difference in the entanglement spec- 
trum can be characterized by the Renyi entropy (37l | 
= jzr^ log Trp" , which reduces to the EE (or von Neu- 
mann entropy) at a — >■ 1. According to Eq. ^ the Renyi 
entropy of Kitaev model is given by Sa — Spa + Sca for 
any a, with Soa and Spa the contribution from Z2 gauge 
fields and fermions, respectively. From Tr^^G[/9^Q] = 
2-(L-i)(n-i)^ one can see that Soa = Sq = (L-l)log2. 
Thus the TEE in Sca is ol independent, which is a generic 



property of the string-net models |38l. l39l|. The a depen- 
dence of Spa in the Abelian and non-Abelian phases 
has qualitative difference due to their different entan- 
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FIG. 3: (Color online) (a) Renyi entropy Sa, and (b) capacity 
of entanglement Ce defined by Eq. ((8]), of non-Abelian (red 
solid line) and Abelian (blue dashed line) states. The black 
dotted line is a linear fitting. The parameters are the same 
as those in Fig. 2. 



glement spectra. If we define p = e~^, the quantity 
50.(1 — 1/a) = logTre^"^ is the same as the free en- 
ergy of a thermal system with Hamiltonian Ji and tem- 
perature t = 1/a. The behavior of the low energy spec- 
trum of T-L can be obtained from the following quantity: 



CE{t)=-t^ [{l-t)S,/t\, 



(8) 



which is termed as capacity of entanglement" and is the 
analog of heat capacity Cy in a thermal system. The 
explicit expression of Sa and Csit) is given in the Sup- 
plementary material, which leads to the numerical results 
shown in Fig. [3] As expected, in the limit of t — >■ 0, Csit) 
vanishes exponentially for Abelian phase but linearly for 
non-Abelian phase, since the latter has a gapless entan- 
glement spectrum with constant density of state. More 
generically, if the entanglement Hamiltonian H describes 
a (1 + 1) dimensional conformal field theory (CFT) in 
long wavelength limit 0, [13] j the capacity of entangle- 
ment is given by Csit) = {TTcL/3v)t for i — >■ 0, with L the 
length of the boundary, and c and v the central charge 
and velocity of the CFT, respectively f4d\. Moreover, if 
H describes a critical theory with dynamical exponent z, 
from dimensional analysis one can obtain the asymptotic 
behavior Csit) oc Lt^/^ for t -J> 0. Thus we see that 
the capacity of entanglement characterizes some impor- 
tant qualitative behavior of the entanglement spectrum 
in generic systems. 
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SUPPLEMENTARY MATERIAL 
A: Derivation of Eq. (4) 

To obtain pA^ we trace out the spin degrees of freedom in B. For the spins in B away from the partition boundary, 
the trace can be carried out in the Hilbert space of fermions and gauge field respectively. However, the spins on the 
boundary sites need more careful treatment. Suppose the links across the boundary are denoted 0,s djibji^ ft — - 
as shown in Fig. 1. (Note that we implicitly assumed the boundary length is even here. The odd boundary length 
case need some extra care but our results below remain valid.) The gauge field Uo„6„ = *7q„7"^ residing on a„6„ has 
spin degrees of freedom from both A and B. To trace out B we introduce new Z2 variables WA,n = *7a2„-i^a2n ^^"^ 
WB,n = il}^^ J ' which are defined on the dotted links in Fig. 1 and belong to either A 01 B exclusively. The 
eigenstates of WA,n, wb.u are related to those of Ua2n-ib2n-n ''^a2nb2n by a unitary transformation. Since \u) is a direct 
product of gauge fields, we denote \u) — \uA,UB,Up) , where |ma); I^b), and \up) are gauge fields in A, B, and on the 
links across the partition boundary respectively. In term of eigenstates of wa and wb defined above, we obtain 

\up) = \wa,wb), (A1) 

wa=wb = {±1} 

if all boundary links have eigenvalues Ua„6„ = 1- For more generic values of Ua^b^, the only change are the signs of 
the terms in the right hand side of Eq. (jAip . which does not enter the reduced density matrix we are interested in. 
Here waib) — {wA(B),n\ and wa = wb means wa.u = WB,n for all n — 1, L. 
In the new basis, the physical eigenstate is given by 

1^) - /^IvTIT f 51 Dg\uA,WA-:UB,WB)\(t'{u)) , (A2) 

^ g,WA=WB 

from which one can obtain the reduced density matrix pA and compute Tr [p^] . 

We define Xg = ilsKlsl-i)/^ J^^g^ 7j tJtI ^'^d = i\9\{\g\-^)/'i Iljeg Vj^ where \g\ denotes to the number of sites in g 
and the ordering of sites in the two products is implicitly taken to be the same such that XgYg = Dg. Xg and Yg are 
the gauge transformation operators that only act on the gauge fields and Majorana fermions, respectively. We further 
denote Xg = Xg.^gnBXg^^gnA and Yg = i\o-\\9B\Yg,Yg,, where Xg, = ^\^-\^\^'^\~^^'^Y{^eaA ^^^^t ^nd 

Ygj^ — il\9A\{\9A\-^-)|^2■ J|.g^ ^ are the gauge transformation operator acting on A region, and similar for Xg^ and Yg^ . 
The ground state (5) is written as 

" /gjV+L+l XI ^g\uA,WA;UB,WB) ■Yg\<i)(u)) , 

^ g,WA=WB 

= /^N+L+i ^9B \UB, W) ■ Xg, \UA, W) ■ Yg,Yg^ \(j){u)) . (A3) 

g.w 

In the last line we have denoted wa = wb as w. The reduced density matrix pA is expressed as 
PA = TrB[|*)(^'|] 

= ^ivTi+T E k^B.w'\xl,Xg, \uBMXgA {ua.w'\xI,^-Ttb,f \YgAYgB \'i^{u)) (</.(«)! yj^yy .(A4) 

g.g' ^w^w' 

To make the interproduct {ub, w'\ A^, Xg^ \ub, w) nonzero, the gauge transformation Xg^ and Xg'^ must be either 
identical in B region, or different by a gauge transformation on all sites in the B region. We define Xb as Xg with 
g = B, and similarly for Xa and Ya.b- Note that Aj^ = (— l)l^^lAg^. So it is consistent to define B — gB such that 
Xg_g^Xgg = Xb for any gB, where Xb = Xg^B- Consequently, we have 

{uB,w'\X^^,^Xgg \UB,W) ^Snj,w' {Sg'j^,gB + ^siwjSg'^^B-gB) (A5) 

with 

L 

Xb{w) = {uB,w\X^g_g^Xg^\uB,w) = {ub,w\Xb\ub,w) = Uij- Wn. (A6) 
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Using this result pA can be simplified to 

PA = ^IvTITT E Y.^9AWA^w){uAMxl,^-Yg,T^B,Fmu)){m\{^ + yB)]Yl 



9 a' 



2Na 



gA,gA-' 



\<f>iu)) {4>iu)\ 



1 + XBiw)YB 



9a 



(A7) 



Notice that Yb — il^KI-^l ^■'^^ Hjes '7j actually the ferinion number parity operator in B region, we see that the 



operator 



1+xb(w)Yb 



is actually the projector of the fermion number parity Yb ~ xb{w). For convenience, we define 



p; 



Xb{w) _ 1 + xb{w)Yb 



B.F 



and similar for A region. The reduced density matrix 



Tr_B,F 



xb{w) 

Pa.f 



(A8) 



(A9) 



is the fermion density matrix with restriction to the fixed fermion number parity xb{w) in the B region. By using 
similar interproduct formula (|A5P for A region, p\ can be computed as 



2Na+ 



1 \ ' -l^ I \ / I V Xb(w) j-tXAit") Xb(w)~,^ 

V^lL^ ^9A \UA, W) {UA, w\ Xjj,^ ■ YgAPA.F Pa,F Pa,F ^g'A 



(AlO) 



gA.gA."^ 

The same calculation can be repeated to obtain 



Pa 



2Na 



Finally we take the trace of Pa by using the same Eq. (jAsp to obtain 



(A12) 



For Eq. (jA6p . we see that flipping the sign of one link Wi, for some integer 1 < i < L, changes the sign of both 
XAiw) and xb{w)- If we define ni7eA(B) ~ Pa{b)i then xa{b){'w) — Pa(b) nn=i ""^"j ^^'^ summation over all 
2^ configurations of w = {wi, i — 1,2, L} leads to 



1 



2(„-l)(L-l) 



Tr, 



( pPA PB V I / p-PA -P 

Ya.fPa.f J + \^A,F Pa, 



(A13) 



In the main text we have discussed that the ground state \u) \4>{u)) must satisfy the constraint D \u) \4'iu)) = 
\u) \(j){u)). Since D = XaXbYaYb and XaXb \u) = paPb \u), we obtain YaYb \4'{'^)) — PaPb \4'i'U'))- In other words, 
the total fermion parity is fixed, so that Pa\Pa^f I't'i'u-)) = Pa^fPa^f — 0. Thus we have 



p-PA _ n 

^a,f — ^ 



(A14) 



so that 



1 



rTr, 



2(n-l)(L-l) 



{pZfPXf) + {Pa'f' pTf) 



(A15) 



where pa,f — p\ f ^ Pa f ^ '^''s [|'/'(''^)) {4'{f^)\] is the free fermion density matrix without fermion number parity 
constraint in the B region. 

To understand this result more intuitively, we note that 



T^^'a,g[Pa,g] = 



1 

2(«-l)(L-l) 



(A16) 
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where pa,g = TrB[|G(M)) (G(m)|] is the reduced density matrix of a pure Z2 gauge field, and the ground state of the 
Z2 gauge field \G{u)) is given by a equal weight superposition of all the 2^"^^ gauge field configurations \u) that are 
gauge equivalent to \u), i.e., \G{u)) — 2^'^~^''/^ \^)- follows that Eq. (jAlsp can be written as 

TtaIpVi - -TtaApXg] ■ "^^aApXfV (A17) 
Thus, Eq. (4) in the main text is proved. 



B: Renyi entropy and entanglement capacitance of Majorana fermions 

As shown in Rcf. [3f], the free fermion reduced density matrix always have the form 



exp 



(Bl) 



with 57 — Y\^^ (1 + 6^*^") the normalization constant and 7„ fermion annihilation operators in the diagonal basis of 
the density matrix. For a system of Majorana fermions 77^, define the equal time correlation function in the ground 
state Cij = {rjirjj) /2 for i^j restricted to the region A. The EE between two parts A and B is then given by 



-Tr[ClogC+(l-C) log(l-C)] 



(B2) 



the 



with C the matrix with entries of Cij. The "single particle density matrix" C plays the role of (l + e' 
thermal entropy, with h the single particle Hamiltonian and the inverse temperature /3 = 1/T = 1. The factor 1/2 
comes from the fact that a Majorana fermion has half the degree of freedom of a Dirac fermion. 

The eigenvalues A„ of the single particle density matrix C in Eq. (IB2[) and e„ in Eq. (jBl[) are related as 



An — 



f 



A. 



e^" + 1 1 - A„ 

Note that for Majorana fermion only e„ > states are summed over in Eq. (jBl|) . Using Eq. (jBip we have 



(B3) 



(B4) 



£„>0 



so that 



Sr. 



1 1 

2 1 -a 
1 1 



^[log(l + e-"^")-alog(f + 6--^")] 

n 



(B5) 



The "capacity of entanglement" Ce defined in Eq. (11) can be obtained by 

t &^ 



2 9t2 

-y 

+2 



t^log Ay* + (1-A„) 



1/* 



2t2 



2cosh(e„/2t) 



(B6) 



It should be noticed that all the results above are for free Majorana fermions. For free complex fermions the only 
difference is an additional factor of 2. 

In the torus case we studied, ky is a good quantum number, and we have 



Ce 



2t' 



E 



2ti V2cosh(e„(fcj,)/2t) 



(B7) 
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For the Abclian phase e„ has a gap since A„(fcj,) does not cross 1/2. If the gap is Eg = min(|e„(fcj,)|), the asymptotic 
behavior of Ce at t — >■ is C'e — e~^o/*' /t^. For the non-Abehan phase, as shown in Fig. 2(b) there is a gapless branch 
of \n{ky) crossing 1/2, which corresponds to ^n{ky) crossing 0. Near ky = Q the asymptotic behavior of |e„(fcy)| is 
|e„(fc,y)| ~ w|fcj,|, which gives Ce oc t/v. In the same way as the heat capacity in a thermodynamic system, CeI^ is 
constant at t — >■ hmit and is proportional to the density of state 1/v at "maximally entangled limit" e„ — >■ or 
A„ ^ 1/2. 



